6. ARIMA Modeling
(6-2) B-J DIFFERENCING
Nao Mimoto - Dept. of Statistics : The University of Akron
TS Class Web Page – R resource page
1. Box-Jenkins Differencing Method
This is called ARIMA modeling.
a) Backward and Difference Operator
Recall the backward operator \(B\), \[ B \, X_t = X_{t-1}. \]
Then define difference operator \(\bigtriangledown\) (del), \[ \bigtriangledown = (1-B). \]
For example, \[ \bigtriangledown X_t = (1-B)X_t = X_t - X_{t-1} \\ \\ \bigtriangledown^2 X_t = (1-B)(1-B) X_t = X_t - 2X_{t-1} + X_{t-2} \]
b) Differencing Linear Trend
Suppose your time series have a line plus a zero-mean staionary noise.
\[ Y_t = a + bt + X_t \] That means \[ Y_{t-1} = a + b(t-1) + X_{t-1} \\ \\ \bigtriangledown Y_t = Y_t - Y_{t-1} = b + X_t - X_{t-1} \]If \(X_t\) is stationary, then \(X_t-X_{t-1}\) is also stationary. Thus we have now \[ \bigtriangledown Y_t = \mu + M_t \] We can try to model this with ARMA(p,q) with intercept \(\mu = b\).
set.seed(42352)
t <- 1:100
Y <- 3 - 10*t + arima.sim(n=100, list(ar = c(.7,.2)), sd=100)
Y2 <- diff(Y) # take the difference
layout(matrix(1:2,1,2))
plot(Y,type="o") # Y is simulated data
plot(Y2, type="o" )
## [1] -9.461098c) Differencing Quadratic Trend
If \(m_t\) is quadratic, then \[\begin{align} Y_t &= a + bt + c t^2 + X_t \\\\ Y_{t-1} &= a + b(t-1) + c (t-1)^2+ X_{t-1} \\\\ Y_{t-2} &= a + b(t-2) + c (t-2)^2 + X_{t-2}. \end{align}\]
That means
\[\begin{align} \bigtriangledown^2 Y_t &= Y_t - 2Y_{t-1} + Y_{t-2} \\\\ &= 2c + X_t - 2X_{t-1} + X_{t-2} \end{align}\]
If \(X_t\) is stationary, so is \(X_t - 2X_{t-1} + X_{t-2}\).
If \(m_t\) is polynomial of deg \(k\), with coefficients $c_0, c_1, $ then applying \(k\) power of difference operator will remove the trend. \[ \bigtriangledown^k Y_t = k! c_k + \bigtriangledown^k X_t. \]
Then you will end up with some statonary series \(\bigtriangledown^k X_t\) with constant trend \(k! c_k\).
t <- 1:100
tsq <- t^2
Y <- 3 - .5*t + .1*tsq + arima.sim(n=100, list(ar = c(.7,.2)))*10
layout(matrix(1:2, 1, 2))
plot(Y,type="o")
plot(diff(Y),type="o")
Box-Jenkins Differencing Turns many non-stationary trends into stationary.
If you difference linar trend \(a+bt\) once, it becomes constant trend \(b\).
If you difference quadratic trend \(a+bt+ct^2\) twice, it becomes constant trend \(2c\).
2. Difference on Random Trends
a) Random Walk
Suppose we have a model with random trend \[ Y_t = M_t + X_t, \] where \(X_t\) is a stationary series and \(M_t\) is a random walk: \[ M_t = \sum_{i=1}^t e_i \mbox{ where } e_i \sim_{iid} N(0, 1). \]
With iid errors \((e_1,e_2,e_3,\ldots)\) random walk is generated as
\[\begin{align} M_1 &= e_1 \\ M_2 &= e_1 + e_2 \\ M_3 &= e_1 + e_2 + e_3 \\ M_4 &= e_1 + e_2 + e_3 + e_4 \\ &\vdots \end{align}\]
We can use cumsum() function in R to do this easily.
Mean of RW is \[ E[M_t] = E\Big(\sum_{i=1}^t e_i\Big) = \sum_{i=1}^t E(e_i) = 0 \]
Variance of RW is \[ V[M_t] = V\Big(\sum_{i=1}^t e_i\Big) = \sum_{i=1}^t V(e_i) = \sigma^2 t \]
Example: Simulate RW
#- repeat this 50 times
plot(M,type="l", ylim=c(-40,40), main="50 random walks")
for (i in 1:50) {
e = rnorm(100) #- copy and paste these 3 lines many times
M = cumsum(e)
lines(M)
}
- Note that mean is 0, but variance increase as time goes on.
b) Differencing RW
Since \[ M_t = \sum_{i=1}^t e_i, \] we have \[ \bigtriangledown M_t = M_t - M_{t-1} = \sum_{i=1}^t e_i - \sum_{i=1}^{n-1} e_i = e_t \]
And we know that \(e_t \sim N(0,\sigma)\).
So if \(M_t\) is Random Walk, then \(\bigtriangledown M_t\) should look like iid \(N(0,\sigma)\) noise.
c) RW with Drift
Random Walk with drift \(\delta\) is \[ M_t = \sum_{i=1}^t e_i \mbox{ where } e_i \sim_{iid} N(\delta, 1). \]
Now each step have little more chance of being positive than negative:
If \(M_t\) is Random Walk with drift, then \(\bigtriangledown M_t\) should look like iid \(N(\delta,\sigma)\) noise.
Mean of RW is now \[ E[M_t] = E\Big(\sum_{i=1}^t e_i\Big) = \sum_{i=1}^t E(e_i) = \delta \]
Variance of RW is unchanged: \[ V[M_t] = V\Big(\sum_{i=1}^t e_i\Big) = \sum_{i=1}^t V(e_i) = \sigma^2 t \]
layout(matrix(1:2, 1, 2))
#- RW with drift (delta = .3)
e = rnorm(100, .3, 1) #- 100 random number from N(1,1)
M = cumsum(e)
#- repeat this 50 times
plot(M,type="l", ylim=c(-40,40), main="50 RW with drift .3")
for (i in 1:50) {
e = rnorm(100, .3, 1)
M = cumsum(e)
lines(M)
}
#- RW without drift
e = rnorm(100) #- 100 random number from N(0,1)
M = cumsum(e)
#- repeat this 50 times
plot(M,type="l", ylim=c(-40,40), main="50 RW w/o drift")
for (i in 1:50) {
e = rnorm(100) #- copy and paste these 3 lines many times
M = cumsum(e)
lines(M)
}
d) Testing RW for Drift
Does this RW has drift or not?
Test if mean of the difference is 0.
## [1] 0.3292302 1.0239543 0.2006950- 95% CI for the mean is \(.33 \pm .20\), so the sample mean is significantly different from 0. We conclude that the above RW has a drift.
Random Walk as Trend
- If \(Y_t = M_t + X_t\), then \[ Y_t = \sum_{i=1}^t e_i + X_t \\Y_{t-1} = \sum_{i=1}^{t-1} e_i + X_{t-1} \] That means \[ \bigtriangledown Y_t = Y_t - Y_{t-1} = e_t + X_t - X_{t-1} \] Since \(X_t, e_t\) are stationary, so is \(e_t+X_t-X_{t-1}\).
If you difference Random Walk without drift once, it becomes N(\(0,\sigma^2\)) iid noise.
If you difference Random Walk with drift \(\delta\) once, it becomes N(\(\delta,\sigma^2\)) iid noise.
Ex1: AAPL
Daily adjusted close price of Apple, Inc. Between 10/17/2007 to 8/4/2008 from Yahoo! finance website
library(quantmod) # Do install.packages("quantmod") if not installed
getSymbols("AAPL") #- download from Yahoo!## [1] "AAPL"
## Series: X
## ARIMA(0,1,0)
##
## sigma^2 estimated as 0.3226: log likelihood=-170.65
## AIC=343.3 AICc=343.32 BIC=346.6 # Arima(0,1,0) w/o drift is suggested
# force drift term. Is it significant?
Arima(X, order=c(0,1,0), include.drift=TRUE) ## Series: X
## ARIMA(0,1,0) with drift
##
## Coefficients:
## drift
## -0.0121
## s.e. 0.0402
##
## sigma^2 estimated as 0.3241: log likelihood=-170.61
## AIC=345.21 AICc=345.27 BIC=351.813. Back to Modeling Lake HURON 4 different ways
Initial Set-up
library(forecast)
source('https://nmimoto.github.io/R/TS-00.txt')
D <- read.csv("https://nmimoto.github.io/datasets/lake.csv")
Lake <- ts(D, start=1875, freq=1)
plot(Lake, type="o")
## Warning in kpss.test(A): p-value smaller than printed p-value## KPSS ADF PP
## p-val: 0.01 0.24 0.025- Stationarity tests results are mixed.
a) Direct ARMA fit
Do not difference, force \(d=0\), and look for best ARMA(\(p,q\)).
# find best ARMA(p,q) by AICc
Fit1 <- auto.arima(Lake, d=0, stepwise=FALSE, approximation=FALSE)
Fit1 ## Series: Lake
## ARIMA(1,0,1) with non-zero mean
##
## Coefficients:
## ar1 ma1 mean
## 0.7665 0.3393 9.1290
## s.e. 0.0773 0.1123 0.3861
##
## sigma^2 estimated as 0.4784: log likelihood=-101.09
## AIC=210.18 AICc=210.62 BIC=220.48
## B-L test H0: the sereis is uncorrelated
## M-L test H0: the square of the sereis is uncorrelated
## J-B test H0: the sereis came from Normal distribution
## SD : Standard Deviation of the series## BL15 BL20 BL25 ML15 ML20 JB SD
## [1,] 0.963 0.952 0.934 0.567 0.641 0.89 0.684## Series: Lake
## ARIMA(1,0,1) with non-zero mean
##
## Coefficients:
## ar1 ma1 mean
## 0.7665 0.3393 9.1290
## s.e. 0.0773 0.1123 0.3861
##
## sigma^2 estimated as 0.4784: log likelihood=-101.09
## AIC=210.18 AICc=210.62 BIC=220.48Analysis 1 (direct fit)
auto.arima() chooses AR(2) with min AICc.
AR(2) with constant mean was fit directly to data. With observation \(Y_t\),
\[{\large Y_t \hspace{3mm} = \hspace{3mm} \mu + X_t }\] \[{\large \hspace{20mm} X_t \hspace{3mm} = \hspace{3mm} \phi_1 X_{t-1} + e_t + \theta_1 e_{t-1} }\] \[{\large e_t \sim WN(0, \sigma^2) }\]
b) ARMA with Linear Trend
If we choose that the level is going down linearly, we can try to fit a line with ARMA errors.
Force \(d=0\) (don’t difference)
Force linear regression by xreg=time(Lake)
# fit the residuals with ARMA(p,q)
Fit2 <- auto.arima(Lake, d=0, xreg=time(Lake), stepwise=FALSE, approximation=FALSE)
Fit2## Series: Lake
## Regression with ARIMA(1,0,1) errors
##
## Coefficients:
## ar1 ma1 intercept xreg
## 0.6682 0.3817 55.5443 -0.0242
## s.e. 0.0936 0.1136 18.0324 0.0094
##
## sigma^2 estimated as 0.4612: log likelihood=-98.66
## AIC=207.33 AICc=207.98 BIC=220.2
## B-L test H0: the sereis is uncorrelated
## M-L test H0: the square of the sereis is uncorrelated
## J-B test H0: the sereis came from Normal distribution
## SD : Standard Deviation of the series## BL15 BL20 BL25 ML15 ML20 JB SD
## [1,] 0.985 0.974 0.969 0.227 0.217 0.783 0.668## Series: Lake
## Regression with ARIMA(1,0,1) errors
##
## Coefficients:
## ar1 ma1 intercept xreg
## 0.6682 0.3817 55.5443 -0.0242
## s.e. 0.0936 0.1136 18.0324 0.0094
##
## sigma^2 estimated as 0.4612: log likelihood=-98.66
## AIC=207.33 AICc=207.98 BIC=220.2Analysis 2 (linear trend)
Linear trend was fit to \(Y_t\), then the residuals were fitted with ARMA(1,1).
We observe $ Y_t = a+b t + X_t $ \(X_t\) is ARMA
In other words,
\[{\large Y_t \hspace{3mm} = \hspace{3mm} a + bt + X_t }\]
\[{\large \hspace{20mm} X_t \hspace{3mm} = \hspace{3mm} \phi_1 X_{t-1} + e_t + \theta_1 e_{t-1} }\]
\[{\large e_t \sim WN(0, \sigma^2) }\]
c) Arima(p,1,q)
- Now difference once by setting \(d=1\).
# Fit the difference of Lake with ARMA
Fit3 <- auto.arima(Lake, d=1, stepwise=FALSE, approximation=FALSE)
Fit3## Series: Lake
## ARIMA(1,1,2)
##
## Coefficients:
## ar1 ma1 ma2
## 0.6385 -0.5349 -0.3514
## s.e. 0.1345 0.1445 0.1055
##
## sigma^2 estimated as 0.4812: log likelihood=-99.88
## AIC=207.76 AICc=208.2 BIC=218.02
## B-L test H0: the sereis is uncorrelated
## M-L test H0: the square of the sereis is uncorrelated
## J-B test H0: the sereis came from Normal distribution
## SD : Standard Deviation of the series## BL15 BL20 BL25 ML15 ML20 JB SD
## [1,] 0.979 0.956 0.942 0.524 0.495 0.577 0.677 # same model is suggested without specifying d=1, because
# auto.arima uses KPSS test to decide on d.
auto.arima(Lake, stepwise=FALSE, approximation=FALSE) ## Series: Lake
## ARIMA(1,1,2)
##
## Coefficients:
## ar1 ma1 ma2
## 0.6385 -0.5349 -0.3514
## s.e. 0.1345 0.1445 0.1055
##
## sigma^2 estimated as 0.4812: log likelihood=-99.88
## AIC=207.76 AICc=208.2 BIC=218.02## Series: Lake
## ARIMA(1,1,2) with drift
##
## Coefficients:
## ar1 ma1 ma2 drift
## 0.6911 -0.6271 -0.3729 -0.0241
## s.e. 0.0956 0.1251 0.1157 0.0100
##
## sigma^2 estimated as 0.4663: log likelihood=-99
## AIC=208.01 AICc=208.67 BIC=220.83
## B-L test H0: the sereis is uncorrelated
## M-L test H0: the square of the sereis is uncorrelated
## J-B test H0: the sereis came from Normal distribution
## SD : Standard Deviation of the series## BL15 BL20 BL25 ML15 ML20 JB SD
## [1,] 0.989 0.979 0.973 0.276 0.239 0.465 0.666auto.arima() INCLUDES drift term by default.
Arima() DOES NOT include drift term by default when d=1.
To force drift term in Arima(), use include.drift=TRUE
Analysis 3 (take difference)
Difference of \(Y_t\) was taken. The difference seems to be ARMA(1,2).
In other words, $ Y_t $ is ARIMA(1,1,2), meaning \[{\large \bigtriangledown Y_t = Y_t - Y_{t-1} = X_t }\] \[{\large \hspace{20mm} X_t \hspace{3mm} = \hspace{3mm} \phi_1 X_{t-1} + e_t + \theta_1 e_{t-1} + \theta_2 e_{t-2} }\] \[{\large e_t \sim WN(0,\sigma^2) }\]
d) Arima(0,1,0)
- difference once, and
## Series: Lake
## ARIMA(0,1,0)
##
## sigma^2 estimated as 0.5383: log likelihood=-106.49
## AIC=214.98 AICc=215.02 BIC=217.54
## B-L test H0: the sereis is uncorrelated
## M-L test H0: the square of the sereis is uncorrelated
## J-B test H0: the sereis came from Normal distribution
## SD : Standard Deviation of the series## BL15 BL20 BL25 ML15 ML20 JB SD
## [1,] 0.169 0.167 0.125 0.161 0.098 0.369 0.733Analysis 4 (take difference)
Difference of \(Y_t\) was taken. The difference seems to be WN.
In other words, $ Y_t $ is ARIMA(0,1,0), meaning \[{\large \bigtriangledown Y_t = Y_t - Y_{t-1} = X_t }\] \[{\large \hspace{20mm} X_t \hspace{3mm} = \hspace{3mm} e_t }\] \[{\large e_t \sim WN(0,\sigma^2) }\]
Conclusions
a) ARIMA(1,0,1) with mean is fitting adequately. This suggests that Lake is stationary (at constant mean) without any transformation.
\[{\large Y_t \hspace{3mm} = \hspace{3mm} \mu + X_t, \hspace{10mm} X_t \mbox{ is ARMA(1,1)} }\]
b) ARIMA(1,0,1) with linear trend is fitting adequately. Intercept and slope are significant.
\[{\large Y_t \hspace{3mm} = \hspace{3mm} a + bt + X_t, \hspace{10mm} X_t \mbox{ is ARMA(1,1)} }\]
c) ARIMA(1,1,2) with drift is fitting adequately. Drift is significant, and same as the slope in (b).
\[{\large \bigtriangledown Y_t = Y_t - Y_{t-1} = b + X_t, \hspace{10mm} X_t \mbox{ is ARMA(1,2)} }\]
d) ARIMA(0,1,0) is fitting adequately. This suggests that Lake is RW w/o drift.
\[{\large \bigtriangledown Y_t = Y_t - Y_{t-1} = X_t = e_t \hspace{10mm} e_t \mbox{ is ARMA(0,0) = WN} }\]
(b) and (c) are concluding the same model.
(a) is a weaker version of (b), but intercept and slope in (b) are significant.
(d) is a weaker version of (c), but ARMA parameters in (c) are significant.
The model (b) (which is essentially same as (c)) is the best candidate.
Summary